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ABSTRACT 


Describes a method for pps sampling developed by H. O. Hartley 
(1966). To facilitate drawing samples using this method, a computer 
subroutine is provided that identifies the units in the sample, calculates 
their sampling probability, and summarizes the distribution of this 
probability over the entire population for use in subsequent estimates 
of the population total and its variance. An estimate of the error 
variance is derived in the case that the ratio of measured value to size 
is a polynomial function of size. This error estimate is a multivariate 
analog of Hartley's linear model for the superpopulation. To permit 
sampling without replacement, excessively large units in the population 
are assigned to a stratum to be measured with absolute certainty. 


THE SAMPLING METHOD 


Sampling with unequal probabilities of selecting a sample unit frequently can 
result in increased efficiency in estimating population totals. The technique 
requires that each unit in the population be characterized by a "size" variable that 
is correlated with the variable of real interest to be measured for that unit. If a 
list of the size variable can be obtained cheaply for all units in the population, 
then sampling with probability proportional to size (pps) according to the method 
developed by Hartley (1966) can be more efficient than pps sampling by other methods. 
How much more efficient depends on the character of the functional relation between 
size and the variable of real interest. 


The distinguishing feature of Hartley's method is that the units of the population 
are sampled at uniform intervals of accumulated size from a random start in the list 
of units sorted into order by size. Hence, the sample that is drawn is more likely 
to contain both small and large units than would a pps sample drawn from a randomly 
ordered list. Hartley has provided an expression for the true variance of estimates 
based on this method of selecting the sample. However, the sample estimate of this 
variance that he suggests is based on the assumption that the ratio of the variable 
of interest to the size variable is a linear function of the size variable. Ina 
subsequent part of this paper, the assumption concerning the linearity of this relation 
is relaxed to include polynomial functions of the size variable. 


PPSORT -- A COMPUTER SUBROUTINE FOR 
SELECTING THE SAMPLE 


PPSORT is a FORTRAN IV computer subroutine that implements Hartley's procedure 
for drawing samples from a population list sorted by size. The program is complete 
with its own composite random number generator (Marsaglia and Bray 1968; Grosenbaugh 
1969) and sorting algorithm. The execution speed of the version of PPSORT listed in 
appendix I was markedly improved over earlier versions by adopting a sorting method 
developed by D. L. Shell, as it was programed by Robert M. Russell of the Pacific 
Southwest Forest and Range Experiment Station. 


The pps selection method used by PPSORT is without replacement. Accordingly, 
the number of sample units must be less than the measure of size accumulated over the 
entire population divided by the size of the largest element in the population. The 
logic of PPSORT is such that if this condition is not satisfied for the number of 
sample units requested, then the population is divided into two strata. The largest 
unit in the population is assigned to one stratum in which every unit will be desig- 
nated for measurement with absolute certainty. The other stratum is comprised of the 
remaining units of the population. The number of sample units to be drawn is reduced 
by one, and the inequality is tested again. If the number of sample units is still 
too large, the process is repeated until the inequality is satisfied. Fortunately, 
the expression for estimating the population total can ignore the separation into 


sampled and completely measured strata. If a stratum for the largest units is required 
to satisfy the above inequality, all the units assigned to this stratum can be given a 
sampling probability of exactly unity. Thus, the estimated population total can be 
seen to be the sum of the estimated total for all the sampled stratum plus the total of 
the stratum of excessively large units. 


To provide the proper elements for the variance estimate based on the units in the 
sampled stratum, PPSORT returns values of the total size of the sampled stratum and the 
number of units selected from it for measurement along with their identification and 
probability of inclusion. Certain other values depending on the distribution of the 
inclusion probabilities over the entire population are also saved by PPSORT for later 
use in calculating the variance estimate. 


CALLING PPSORT 


To use PPSORT, a computer program would be prepared that enters the data describing 
the entire population, calls PPSORT, and then displays the identification and sampling 
probabilities for each member of the sample drawn by PPSORT, and the values required 
for the subsequent variance estimate. A simple version of such a program with test 
data that can be used to verify the functioning of the PPSORT program is provided in 
appendix II. Three random odd integers must be assigned to the variables KR, LR, and 
MR to start the pseudorandom number generator used in PPSORT (Marsaglia and Bray 1968). 


The subroutine is executed by the statement 
CALL PPSORT (NP,ID,CHAR,NSAM, IDSAM,PP,TX,TS,TP,NRS,N,LR,MR, KR) 


in which the first four and the last three arguments of the calling sequence are estab- 
lished prior to the execution of the CALL statement. 


NP contains the number of units in the population. (If NP exceeds 498, 
the DIMENSION statement in the PPSORT program should be changed 
accordingly.) 


ID is a vector of NP elements identifying the units of the population. 


CHAR is a vector of NP elements containing the sizes of the units of the 
population corresponding to the same sequence as the ID vector. 


NSAM contains the number of units to be measured (sample units plus those 
that may be designated for measurement with absolute certainty). 


The intervening seven arguments contain information about the sampling as it was 
accomplished by PPSORT. 


IDSAM is a vector of NSAM elements containing the values copied from the ID 
vector for only those units selected for measurement. 


PP is a vector of NSAM elements containing the value of the probability 
with which each unit was selected for measurement. Their sequence is 
the same as the sequence of ID's in IDSAM. (PP(I) = Ps) 


TX contains the sum of the size variable for the stratum to be sampled. 
TS,TP contain components of the sample variance estimate defined on page 4. 


NRS contains the number of units drawn from the stratum to be sampled 


(NRS = 7). 
N contains the total number of elements in the sampled stratum. 
KR 
LR contain odd random integers used by the random number generator. 
MR 


CALCULATING ESTIMATES 


Population Total 


The population total is estimated from a pps sample by 


; Oe ge ee 1 
el ay sae ae 
Z=1 e=1 - 
in which 
n is the number of units measured in the sample. 
c is the number of units assigned to the stratum of too-large 


units to be measured with absolute certainty. 


Y, is the variable of interest for the 7th measured unit. 
NW 
P = 7 x / zr ee the probability with which the zth unit was selected for 
j=l 


measurement where x. equals the measure of size of the 7th unit in 
the sampled stratum. P; equals unity for units in the too-large 
stratum. 
The last summation would not contribute to the variance of J. 
Note that P; depends on m. Accordingly, inadvertent change in after the P. 


is calculated should be avoided. If u sample units are not measured for reasons 
independent of their sizes, then Y could be estimated by 


t=n+1 


in which Y, = 0 for the omitted sample units. 


Variance of Population Total 


The variance ie the population total is estimated from a sample drawn with PPSORT 
by the expression’ 


Var Y = n2 (b2 - Var b) TS + n TP G2 


in which 
n n n 
ina ee WG eR ag P./n 
ps t=1 t=1 g=l1 
n n 
r ae (C2 )2/n 
i=1 j=l 
n n n n n 
Eo AY. /Pe ae PP een = bel wee eh eel 
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n-2 


‘ i n n 
Var b = o2/| = P.2 - @2z Poin 
= r=] 


Sy = probability with which the sample cluster of m units was selected 
for measurement 
N = total number of units in the population to be sampled 
N n N n 
TS =. Ys Ch eh fnes see me 
aie eee gain eye 
N 
TE a lem ie Oe 
ker * 


The last two values are calculated by PPSORT. 


This expression is appropriate when the trend of Y;/P; is linear with respect to 
P;. If the sample indicates that the trend is curvilinear, the variance will be over- 
estimated unless the formula is modified as described in the following section. 


IThe factor (N-1)/N in Hartley's expression is replaced here by the exact value 
(1 - Is, *) because it is readily available when PPSORT is used. 


Derivation of Variance Estimate for Polynomial Model 


The variance estimate given above assumes that the trend of Y./P. is a linear 
function of P.. If the relation is known to be curvilinear in advdncé, a suitable 
transformation of the size variable should restore the linearity. However, if the 
curvilinear trend is not observed until after the sample has been measured, the 
variance estimate based on the linear model would be too large. In such a case, 
the following analysis shows how a multivariate analog of the linear model could be 
used to estimate the variance of the population total. 


Let the following symbols represent vectors or matrices of the indicated 
dimensions. The subscript k designates the kth sample cluster defined by the sorted 
order and the sampling interval. The range of k is from 1 to J. 


= (¥,,/? 3) of dimension xl 
m = (1) W " nNxl1 
R = (m'r)) i " Nx1 
2s TPN Paster ome m " NXqQ 
x = * (m'xX)) " ” Wxq 
W = (8,557) uM m NxN 
M = (1) " Ll NxN 


In the above vectors and matrices, the scalars are defined as follows: 


Y, is the measured value of the variable of interest for the tth unit in 
the population. 


P is the probability with which Y, is measured. 
Sy is the probability that the kth sample cluster of Y; is measured. 
Tes Sea) 
ee Oif iff 
N is the total number of units in the population. 
n is the number of units in the sample. 
q is the degree of the polynomial function of Pp. 


Note that M is a matrix and m is a vector having unity in all elements. The effect of 
multiplying a vector by M or Mm is to replace each element in the vector with the 
column total. In particular, 


N 
MWM =M_ because Zs, =] 
kK 


With this notation, define three additional matrices: 


Z = X - MWX of dimension Nxq 
B = (Z'Wz)~? Z'WR " a lier sl 
EB =rRy = "ZB Lf WY WVx1 


in which Z'’ is the transpose of Z. 
Now Z'WZ is the analog of 


N eee N 
n2 F Ss, (&, - Def es x, )2 
R= B= 
and (R-MWR)'W(R-MWR) is the analog of 
N wi Oa N 
Bee, 3G r)2 (@ez a)? 
k=1 t=1 


in Hartley's notation. 


Furthermore, because MWM = M, it can be shown that MWZ = 0, the null matrix of 
dimension nxq. It is also true that Z'WE = 0, the null vector of dimension 
qxl. 


The objective is to derive an estimate for the variance of the estimate of the 
population total Y = mir. The variance of Y is given by 


Var Y = n? (R-MWR)' W (R-MWR) 
= mn? (ZB + E - MWZB - MWE)' W (ZB+E - MWZB - MWE) 
= n* (ZB + E - MWE)' W (ZBtE - MWE) 
= mn? (B'Z'WZB + (E - MWE)' W (E - MWE) + 2 B'Z'W (E-MWE)) 


The last product is zero because 
BtZtWE-="B!0) = 0 

and 
B'Z'WMWE = B'(MWZ)' WE = B'O'WE = 0 


In the variance expression, the first matrix product, Z'WZ is the multivariate analog 
of the factor 7S computed in PPSORT. 


To estimate the B'Z'WZB component of Var Y requires first an estimate of B based 
only on the sample. Let b represent the vector of g estimates of the elements of B, 
and let z be the matrix of deviations of the q powers of P; for each of the m units 
in the sample: 


n n ; 
= a Dre 2 Cl a7, 
Z (2p y P/n, Bs 2x P. /n, ab 54 P; /n) 


=i Zl B 
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THOSE INTERESTED IN KEEPING ABREAST OF ANY CHANGES, 
ADDITIONS, OR CORRECTIONS IN THE MATERIAL PRESENTED 
IN THIS REPORT SHOULD FILLIN THEIR NAME, AFFILIATION, 
AND ADDRESS BELOW. THE CARD SHOULD BE DETACHED AND 
RETURNED TO: 


Forestry Sciences Laboratory 

Intermountain Forest & Range 
Experiment Station 

Attention; Albert R. Stage 

1221 So. Main St., P.O. Box 469 

Moscow, Idaho 83843 


PPSORT User: RP-INT-88 


Name 
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City State Zip Code 
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Finally, let r be the ~ element column vector of ratios: 


Then define the estimate of B to be 
b = (z'z)7! zr 


Next, define = as the expected value of 
(b-B) 'Z'WZ(b-B) 
Because Z'WZ is a matrix of constants, 
Meoge (b'Z"WZb) —<B*Z"WZE 
in which ce (*) means "the expected value of" (:). 


Under the assumption that r is drawn from a superpopulation having a mean of ZB 
and variance equal to e (e'e) = o%, 


(b-B) = (z'z)7+ z' (zBte) - B 


B+ (z'z)7} z'e -B 


= (z'z)7} z'e 
Then, 
e (b-B)' Z'WZ (b-B) = e (e'z(z'z)"? Z'WZ(z'z)”> z'e) 
To simplify notation, designate the symmetric matrix between e' and e by C, then 
: n ne 21-1 
e(e'Ce) = me € (e,7) Se cs re € (22,0; .) 


The expected value of e,;? is o* and the expected value of e,;e, is zero if the errors 


are independent for a sufficiently high degree polynomial. J 
Therefore, 
‘ n 
= = i 
y= a C. o* tr(C) 


Because one can commute matrices under the trace operator, 
tr(C) = tr (z(z'z)74 Z'WZ (z'z)-* Z) 


Petey ea) 


tr (Z'WZ (z'z)7!) 


Hence, 2 can be estimated by the residual mean square error times the sum of the 
diagonal elements of Z'WZ(z'z)-1. 


The final component of Var Y to be estimated is (E - MWE)' W (E - MWE) which 
equals E'WE - E'WMWE, or returning to the summation notation: 


N ¢-1 N WM t=1 
s £7 #2, & 2h Bs - 8 ,°E * -2 2 2 6b Es. 
1 f=1 j=l Jd fe] f=1 j=l Jd 


a 


Lt 


Taking expectations, and noting that each element E. is the average of residuals, 


ee ga o2 ul 
e (E-MWE)' W(E-MWE) = — 3s. - —¢zs.2 = (—— () 29*rts#2) 
n E n t n D 
t= sl =i 
Putting all the components together, 
ome 215 4" G2 Z 
Var Y = n2 b'Z'WZB - 62 tr [Z'WZ(z'z) —] + a s.*) 
Lai 


in which 
o* = (r"r - bz'y)/ @ =@ = 1) 


and 


When g = 1, the above matrix equation reduces to the variance estimate in the 
previous section for the linear model. Obviously the order g of the polynomial is 
not known when PPSORT is executed, so the gxq matrix Z'WZ must be calculated 
subsequently by a procedure analogous to the computation of 7S in PPSORT. Such a 
program could be prepared by omitting from PPSORT the selection of random sample 
units and by making SX a q-element vector of powers of the difference CP(J+1) - CP(/J). 


EXAMPLES 


The first example shows the output from sampling the same population described by 
Hartley (1966). 


ID ESiie) SLZE NO. OF UNITS TO BE 
MEASURED = 3 
10 130.00000 
8 109. 00CC9 
6 105.00000 
4 9300000 
2 382.00000 
9 115.00000 
7 109.00009 
5 95.00000 
3 93.00CCO 
1 80.00000 
SAMPLES DRAWN 
NOs IN STRATUM = 10) 
NO. OF SAMPLES = 39 


O-LOLLODOE 04) 
0's:338595380E—G3 ) 
0-38516406E 900) 


— eee 
44 
nx 

ou dt 


ID PROBABILITY 


3 0.2759644E 00 
6 O.S1ISTZIE 00 
9 POR S424 57 EL 00 


From the population of 10 units, three are to be drawn for measurement. The third, 
sixth, and ninth units comprise the sample. Table 1 shows the calculation of the 
estimates of the population total and its variance based on the measurements (Y7) of 
these three sample units. The estimate of the population ,total is the sum of Y7/Pz = 
578. The sample standard error of this value is 7642.61 2 on 2) 25)..5: 


2 This value differs from Hartley's calculated variance of 653.2 only because the 
exact value of (1 - = 8,7) has replaced his (W-1)/N approximation. 


Table 1.--Calculation of sample mean and variance estimates (W=10, n=3) 


Unit De ae Ay Ee 
E L Dit 
3 0.27596 56. 202.924 
6 SSIS 7 62% 198 .990 
9 - 34125 60. |W Ass YAs) 
Sums - 9287834 WAS} 5 Si aa As) 
Means - 509593 S9RS555 192.581 


Sum of Squares 0.2896829 10580.00 111690.465 
Y= aE Sah 

L L 
YF Gas (03 P.)?/n = 0.2896829 - (0.9287834)2 / 3 = .0021367334 


zr Y, - = (Y,,/P;) y P ,/n = 178. - (577.74128) (0.9287834)/ 3 = -0.8655035 


M 


(Y,/P,)? =e Y /P,)?/n = 111690.465 - (577.74128)2/ 3 = 428.80126 


Ss 
i 


-.8655035 / .00213673 = -405.059 


> 
Nh 
| 


[428.80126 - (-405.059)2 (.00213673)] / (3-2) = 78.22111 


b? - Var b = 164072.7935 - 78.2211 / .0021367334 = 127464.9976 


TS = .000385958 
TP = 0.85164 
Var Y = 9 (127464.9976) (.000385958) + 3 (0.85164) (78.22111) 


442.765 + 199.849 = 642.614 


10 


The next example illustrates the way in which PPSORT stratified the excessively 


large elements (added to the same population used previously) into a stratum to be 


measured with absolute certainty. 


_ 
© 


Ke ODOANOMNAWN & 


~~ 


Notice that there are three samples, although four units are to be measured. 
Consequently, the error components are exactly the same as in the previous example, 
although the estimate of the population total is increased by the measured value of 


the eleventh sample unit. 


ESie, SIZE 


80.90000 
82.-000C0 
93.09000 
93.00000 
95200000 
105.0CCCO 
109. 00000 
109.090C0 
115.00C00 
130.000CO 
9399.000C0 


NOS OF UNITS TO. BE 
MEASURED = 4 


SAMPLES DRAWN 


NOs IN STRATUM 
NOs OF SAMPLES 
= Qe 1LOLLOOOE 04) 
TS) o> 0 3859580E-03 } 
= 09.8516406E 00) 


“ol 
~ 
(=) 


3) 


—_ —_ 
+ 
~< 


ID PROBABILITY 


lt O.1LOO0000E O1 
4 022759644E 90 
6 0as1TIs5721E 00 
9 0.-3412457E 90 


1a 


to the stratum for measurement with probability equal to one. 


The last example illustrates a situation in which two units must be assigned 


The remaining three 


sample units happen to include all of the remaining "large'' units in the population. 
Nevertheless, the smaller units could have been a part of the sample had the random 
number been small enough. 


— 
3) 
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ESte, SIZE 


3.00009 
3.00000 
5200000 
5.00000 
109.09000 
109.0)000 
115.00C900 
13C.90000 
999.09000 
1.00000 
1.500CO 


NO. 


~ eR RR 


OF UNITS FG BE 
MEASURED = 5 


NO. 
NO. 
TX 
TS 
TP 


SAMPLES DRAWN 


IN STRATUM = 9) 
CF SAMPLES = 3) 
= 0-3515000E 03) 
7 0.1206964E-01) 
= 0.2865710E OO) 


ID PROBABILITY 


11 O.1OQ00000E Ol 
10 O.1LOOODDOE OL 
t— O.9302982E (00 
8 0.9302988E 00 
9 0.98 150331 00 
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APPENDIX | 
PPSORT Program 


SUBROUTINE PPSIRT(NNNe I DeCHAR, NNSAM,ZIUSAM,PPeTXeTSeTPeNSAMeNygLRy 
1 MRsKR) 
DIMENSICN VEC (500) ,CP(500) ,IND(500),CCP(500) 
THE ABOVE FCUR VECTORS MUST BE CIMENSIONED GREATER THAN NNN + 2 
DIMENSION ID(VNN) ,sCHAR(NNN) ,1TDSAM(NNSAM), PP(NNSAM),NRE128) 
EQUI VALENEE { VEC,.CCP) 
MLReMMReMKR, AND NR ARE USED IN THE RANDOM NUMBER GENERATOR. 
TO START THE RANDOM NUMBER GENERATOR, THREE ODD RANDOM INTEGERS 
SHOULD BE ASSIGNED TO LRe KR» AND MRy RESPECTIVELY. 
SEE DHE REPERENGE CLiRED BELOW. 
DATA MLR»s MMRe MKR / 65539, 33554433, 362436069/ 
DATA NR / 
PPlig, SOlte LOOsme 2 Ollie 4455) 21S 4159 430— 277s 3539 153% 237s §4654 
Olisn 2O ss) 415s 4415 S69, 995 D5Dil> 4635 LTS, 555 “4833's 7954 "3635 
“Big COB es Lie 4095 BOS 2315 43'99 209, 357% 1435 2055 255 Ils 
Sills e486 599s) Silke 4539 Oil's Zils BOTs Sige L235 9243's 29; 
Sa we Zlihis pS Ol, B5i9 S30 229 4 493, 2633 BIT, 153 449, 45% 279, 
BZSie0 ON) LO Ms Lis el Sos) Sis) OOO es) LOis LO5Ss L3i%s “1215 3655 
B3a 289s 41734 4055 SOs FIs 4EBle 3B51y 1595 3495 295, 49, 359 
393% 19+ 265% 285% Ws Vows 45s “479%, LOZ, 25% 4595 935 "3316 
LE9ne ZOD Smo 2s 2s 2596) 335i Lids 2335 145, Lits 139% 9, 431, 
BO ee 4g 4451¢) 295 1S, Tl, SES, 85s B20 Ss 63, 379 
N = NNN 
NSAM = NNSAM 
THIS PROGRAM WILL DRAW 'NSAM* SAMPLES FROM TFE* ID*' POPULATION OF 
*N*® UNITS WITH PROBABILITY PRUPCRTICNAL TC TFE ASSOCIATED VALUE 
OF "CHAR* WITHOUT REPLACEMENT. 
THE SAMPLE IS ORAWN SYSTEMATICALLY FROM THE LIST ORDERED BY 
INCREASING VALUES OF *CHAR*t, USING THE THEORY CESCRIBED IN: 
HARONMEY¢ Heer VOURS (AMER. STAT. ASSOC. 61(215)3 7139-748. 
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THE UNITS SEEECTED WILL BE IDENTIFIED IN THE VECTOR * IDSAM*. 


"NNSAM* MAY BE DECREASED SO THAT THE STEP IS GREATER THAN THE 
MAXIMUM VALUE OF "CHAR'. IN THAT EVENT, 

THE LARGEST ELEMENT OF THE POPULATION WILL BE ASSIGNED A SAMPLING 
PRUBABILITY OF UNITY AND DELCTED FROM THE POPULATION. 


TOP SSe ide oer SUE USED IN, THE SYSTEMAT TEC SELECT LON. 


FOR PURPOSES OF ERROR CALCULATICN, THE SUM OF SQUARED DEVIATIONS 
OF THE MEAN OF CHAR FOR EACH OF THE "N® POSSIBLE SAMPLES 

FRCM THE WEIGHTED MEAN OF THe SAMPLE MEANS, EACH WEIGHTED BY THE 
PROBABILITY WITH WHICH THAT SAMPLE CLUSTER MAY BE SELECTED, IS 
CUMPUTED IN *TS*. ONE MINUS THE SUM CF SQUARES OF THE PROBABILITY 
PSPCALCULAMED. TIN) Sa) P' 

THE REMAINING PARTS OF THE EXPRESSICN FOR SAMPLING ERROR MUST BE 
CALCULATED AFTER THE SAMPLE IS* MEASURED. (SEE HARTLEY*S £Q.(29)). 


IF ( N «LT. NSAM) NSAM = N 

iC UNSAM SEES 10) RETURN 

VEC(1L) = 0.0 

IND(N+1) = 1 

CHAR({N+1) = 0.0 

FIRST SORT INT) INCREASING ORDER CF CHAR WITFKCUT RE-ARRANGING 
BETZ =" Tsitt 

IND(J) = J 

NA = Nk 
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16 


20 


23 
25 


30 


IF (NM1L) 49495 
VEC(2) = CHAR(1) 


CE Gah atitiee aN, GO iG eZ. 
MM) S=UMAXOL My M2) 1) 
I = 


IM=I1+4 
J =o 
JM = IM 


IHOLD = IND(IM) 
JR = IND(J) 
LEC HARGIRY .Sleh. CGHARICLHORD i eGo wir Ss 
IND( JM) = IND(J) 
JM = J 
i=) 4 — 7M 
VEO SG. OMe GOehOr ES 
INO( JM) = IHOLD 
ely tl 
IM = [ *°M 
TE OEM Been) GO Taal 
M = M/2 
CE MaGill On) COs TORS 
DO 20 J = 1 oN 
VEC(J+1) = VEC(J) + ABS(CHAR(CINO(JS))) 
NOTE *VEC* IS INDEXED PLUS ONE) 
NOW COMPUTE TRIAL STEP. 
b= 
IMAX = INDI{N) 
D = VEC({(N+1L)/ FLOAT (NSAM) 
IF ( DsGTs.CHAR(IMAX)) GO TO 30 
IDSAM(T) = IOCIMAX) 
PIPE = le10 
Teale 
NSAM = NSAM-1 
N=N-1 
LE A NSAM, SGT. VO De GOMnien23s 
De= VECO) 
TS = 0. 
BP =~ 16.0 
RETURN 
CONTINUE 
UNIFORM RANDOM NUMBER GENERATCR FOR IBM 36C 
FOR OTHER COMPUTERS,sSEE THE FOLLOWING 
REF: MARSAGLIA AND BRAY, CCMM. ACM 11(11), 757-759 
GROSENBAUGH, Le ResCOMM. ACM 12(11),639 


= LR*MLR 
MR = MR¥*MMR 
= 1 + LPABSTERY. ROT T2106 
SEL=(0.5 + FLOAT(NR(JR) + LR + MR) *0.23283064E—-9) *D 
KR = KR*¥MKR 
NR(JR) = KR 


FINO SAMPLE UNITS. 
FNS = FLOAT(NSAM) 


CPi) = 10210 

DO 40 J = l»N 

CP(J+1) = FNS¥VEC(J#1) /VEC(N#1 ) 

CEPI = Cet(Iel)— BLOATCIEIX(EPLILI)) 
IF «( SEE.GT.VEC(Jtl)) GO TO 40 

FOUND A SAMPLE UNIT. 

JJ = IND(J) 

IDSAM( I) = ID(JJ) 

PPT) — CeCJEt)— Ceti 


SEE-= SEL + D 
NOW CALCULATE THE ERROR COMPONENTS DEPENDING CN DISTRIBUTION OF 
*CHAR® THROUGHOUT THE POPULATION. 
40 CONTINUE 
TX = VEC(N+1) 


TSS = 0.0 

iSs7= 0.0 

SEL = FNS/FLOAT(N) 
SX. => €P (C2) 

SP = 0.0 

TP. = 1.0 

a=. tk 


CCP(N) = 1.0 
CALCULATE SUM SF P(1) FOR FIRST SAMPLE CLUSTER 
DO 120 J. = 2yN 
LE Me leGilie LELXVERK I+ 1)))) } GO 10, 120 
SK sox EE PIE) = CP Us) 
Pe=yl ty 
le { [ <GESNSAM 3) GO TGQ 130 
12C CONTINUE 
130 CONTI NUE 
DO 180 J = 1ly»N 
TRY = 6-E+35 
DO 160 I = 1,N 
le” CEP(T) .GeE. TRY ), GO TD) 160 
TRY = CCP Cl) 
IMIN = I 
IMIN INDEXES UNIT TQ DROP FROM PAST SAMPLE CLUSTER 
160 CONTINUE 
Pa — CEPLIMIN) — SP 
SP = CCP(IMIN) 
CCP(IMIN) = 8.£+35 
XBAR = SX/FNS — SEL 
TSS = TSS + P¥*¥XBAR*¥XBAR 


TS = TS + P¥*XBAR 
TP = TP — p#pP 
180 SX = SX + CP(IMIN) + CPC(IMIN#2) - 2.0*CP(IMIN+1) 
1So=7SS.— 1S40S 
RE TURN 


END 
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APPENDIX II 


Calling Program Example 


DI MENSI GN PP(20) sIDSAM(20) ,1D(200),CHAR( 200) 


KR = —32860S067 

LR = 1409859205 

MR = 402656419 

WRITE (64906) 
906 FORMAT(1H1) 


10 READ (5+908,END=200) NP,NSAM 


908 FORMAT(3I15) 


12 READ (5,909) (ID(1I),CHAR(I) »[=lyNP) 


909 FORMAT(15,F10.5) 
WRITE (63904) NSAM, 


( 


ID(T),+CHAR(1),1=1L,NP) 


904 FORMAT(///14X5*POPULATION'//1H 8X_,'ID*,1LX,*EST. SIZE"6X,*NO.w OF U 
INITS TO BE*/40X, "MEASURED =", 1570111, F20.5))) 
2C CALL PPSORT(NP sIDsCHARsNSAMsIDSAMePPy TX oTS eT Pe NRSoyNoL Ry MRe KR) 
IF ( NSAM «GT. NP) NSAM = NP 
WRITE {65907) NsNRSoTXeTSeTPe( IOSAM(1I),PP(1I),f=1l,NS AM) 
907 FORMAT(46X*SAMPLES DRAWN'//40X*( NC. IN STRATUM =", 15,')*/40X, 
LA NOS OF SAMPLES) 3g iltd) gi OH GOX Ge tl pK eGo lig ex Xie nee Neary 


2 EMG lise Ne FOX SEP. 
ZIM 2Xceul Sil Sent dh) 
READ (5,908,END=200) 


= FEMG 9 8g 45X eo TOM 2X3 UPROBAB ILE, 7, 


NP,»NSAM,NEW 


DEAS ANEW te Es OP GOs an Grak2 


WRITE (64,910) NSAM 
910 FORMAT('O ADDITIONAL 
IATION.'//) 
GO TO 20 
20C RETURN 
END 


",I8,' UNITS TO BE MEASURED FROM SAME POPUL 
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